%% Simulate Nonlinearities During Disinflation
% by Jaromir Benes
%
% Simulate the nonlinear credibility model to show different outcomes
% during disinflation, depending on the initial level of credibility.
% Evaluate the effect of nonlinearities compared to a linearized model.

%% Clear Workspace

clear;
close all;
home;
irisrequired 20140401;
%#ok<*NOPTS>

%% Load Model Object

load read_model m;

%% Simulate Permanent Unanticipated Disinflation
%
% Simulate an unanticipated disinflation shock (reducing the inflation
% target by 1 pp) in four different ways:
%
% # With the first-order approximate solution (linearized model), i.e. when
% neither nonlinear Phillips curve nor the credibility matter.
%
% # With taking into account both nonlinearities, starting with the
% central bank fully credible, $c_t = 1$.
%
% # With taking into account both nonlinearities, starting with low
% credibility, $c_t = 0.5$.
%
% # With taking into account both nonlinearities, starting with low
% credibility, $c_t = 0$.
%
% Change `phi=1` <?phi?> to allow for permanent changes in the target
% through the shock `et`. Note that this parameter change does not affect
% the steady state (no need to call the function `sstate` again) but it
% does change the solution matrices, and hence the model has to be resolved
% <?resolve?>. Create a database with steady-state time series for each
% variable, and create a shock to the target <?sstatedb?>.
%
% Simulate first the linearized model; in this case, the Phillips curve
% convexity and the credibility do not matter <?linear?>.
%
% Run a nonlinear simulation starting with full credibility <?fullCred?>.
% This is the default, steady-state value found in the input database `d`
% -- hence, the initial condition for `c` does not need to be changed
% <?noChange?>. To run a nonlinear simulation, use the option
% `'nonlinear='` <?nonlinearOpt?> assigning it a time horizon over which
% the nonlinearities will be taken into account in each period (beyond that
% horizon, the linearized solution is assumed).
%
% Simulate again, starting with low credibility; adjust the initial
% condition for `c` accordingly <?lowCred?>. Simulate one last time,
% starting with no credibility. Change the initial condition for `c` to
% zero <?noCred?>. In addition, because the nonlinearities kick in more
% severly, increase the maximum number of iterations (`'maxiter=' 5,000`
% instead of the default 100), and reduce the step size in each iteration
% (`'lambda=' 0.5` instead of the default `1`). Otherwise, the simulation
% would crash (try that by changing these options).
%
% Nonlinear simulations return three extra output arguments:
%
% * A flag (`true` or `false`), `flag2`, `flag3`, `flag4`, indicating
% convergence.
%
% * A time series of add-factors, `af2`, `af3`, `af4`, that need to be
% added to the first-order approximation of nonlinear equations (i.e. those
% marked with `=#` in the model file) to achieve an exact nonlinear
% simulation.
%
% * A time series of discrepancies between the LHS and RHS in nonlinear
% equations, `di2`, `di3`, `di4`.
%
% These extra output argument are analyzed below.

m.phi = 0; %?phi?
chksstate(m);
m = solve(m); %?resolve?

% d = sstatedb(m,1:40); %?sstatedb?
% d.et(1) = -1;
% 
% s1 = simulate(m,d,1:40); %?linear?
% s1 = dbextend(d,s1);

% % d.c(0) = 1; %?fullCred? %?noChange?
% [s2,flag2,af2,di2] = simulate(m,d,1:40,'dbOverlay=',true, ...
%     'nonlinear=',40,'tolerance=',1e-8); %?nonlinearOpt?

% d.c(0) = 0.5; %?lowCred?
% [s3,flag3,af3,di3] = simulate(m,d,1:40,'dbOverlay=',true, ...
%     'nonlinear=',40,'tolerance=',1e-8);

% d.c(0) = 0; %?noCred?
% [s4,flag4,af4,di4] = simulate(m,d,1:40,'dbOverlay=',true, ...
%     'nonlinear=',40,'tolerance=',1e-8,'maxIter=',1000);

%%

T = 40;
d = sstatedb(m,1:T);
d.c(0) = 0.5;
d.ey = tseries(1:T,@randn)*0.05;
d.epi = tseries(1:T,@randn)*0.05;
d.pi(0) = 5;

%%

s0 = simulate(m,d,1:T,'anticipate=',true,'nonlinear=',[]);

s1 = simulate(m,d,1:T,'anticipate=',true,'nonlinear=',60,'maxiter=',1000, ...
    'display=',100,'lambda=',0.7748);

%s2 = simulate(m,d,1:T,'anticipate=',true,'nonlinear=',T,'maxiter=',1000, ...
%    'sparseShocks=',true);
return

%%

p = plan(m,1:T);
p = exogenize(p,{'y'},1:T);
p = endogenize(p,{'ey'},1:T);

ss1 = s1;
ss1 = rmfield(ss1,'ey');

z0 = simulate(m,s0,1:T,'anticipate=',true,'nonlinear=',[],'plan=',p);

%%

TT = 5;

[z1,flag,addf,disc] = simulate(m,ss1,1:TT,'anticipate=',true,'nonlinear=',TT, ...
    'maxiter=',100,'plan=',p,'sparseshocks=',false,'lambda=',1, ...
    'ignoreshocks=',false,'display=',1, ...
    'solver=','plain');

return

%% Report Results
%
% Plot the model's five variables toghether with the sacrifice ratio,
% i.e. the annualised cumulative output loss.

dbplot(s1 & s2 & s3 & s4,0:40, ...
    {'"Output" y','"Inflation Q/Q" pi','"Inflation Y/Y" pi4', ...
    '"Policy rate" r','"Credibility" c','"Sacrifice ratio" cumsum(y)/4'}, ...
    'tight=',true,'zeroLine=',true,'subplot=',[3,3]);

grfun.bottomlegend('Linearized', ...
    'Nonlin, Full Init Credibility', ...
    'Nonlin, Low Init Credibility', ...
    'Nonlin, Zero Init Credibility');

grfun.ftitle('Disinflation');

%% Nonlinear Simulation Diagnostics
%
% Check the two diagnostics time series returned from nonlinear
% simulations:
%
% * The add factors in the nonlinear equations should be dying out towards
% the end of the simulation <?plotAddF?>
%
% * If the simulation converges successfully, the discrepancies are lower
% than the tolerance level set by the `'tolerance='` option <?plotDisc?>
% (if not specified, the default tolerance level is `1e-5`).

figure();

subplot(3,2,1);
plot(af2{1}); %?plotAddF?
grid on
title('Add-Factors in Simulation #2');

subplot(3,2,2);
plot(di2{1}); %?plotDisc?
grid on
title('Discrepancies in Nonlin Equations in Simulation #2');
grfun.yaxisreformat(gca());

subplot(3,2,3);
plot(af3{1}); %?plotAddF?
grid on
title('Add-Factors in Simulation #3');

subplot(3,2,4);
plot(di3{1}); %?plotDisc?
grid on
title('Discrepancies in Nonlin Equations in Simulation #3');
grfun.yaxisreformat(gca());

subplot(3,2,5);
plot(af4{1}); %?plotAddF?
grid on
title('Add factors in Simulation #4');

subplot(3,2,6);
plot(di4{1}); %?plotDisc?
grid on
title('Discrepancies in Nonlin Equations in Simulation #4');
grfun.yaxisreformat(gca());

grfun.bottomlegend('Phillips Curve','Credibility Signal Equation');

grfun.ftitle('Nonlinear Simulation Diagnostics');

%% Evaluate Discrepancies between LHS and RHS in Each Equation
%
% Evaluate the discrepancies between the LHS and RHS for each equation.
% This can be, of course, done only for databases simulated in a full level
% mode (not deviations).
%
% The function `lhsmrhs` evaluates the discrepancy between the LHS and RHS
% in equations by substituting the actual future observations for the
% variables, their lags and leads. The results are therefore only valid
% when all shocks are anticipated. The output argument `D` is a NEq-by-NPer
% array, where NEq is the total number of equations in the model and NPer
% is the number of periods where the discrepancies are evaluated. The order
% of equations in the individual rows of `D` is the same as they appear in
% the model file, `credibility.model`.

D = lhsmrhs(m,s3,1:40); %?evaluateLhsmRhs?
size(D)
D(:,1:5)

figure();
subplot(1,3,1:2);
plot(D.');
axis tight;
grid on;
grfun.yaxisreformat(gca(),'%.1e');
title('Discrepancies between LHS and RHS in Each Equation');

labels = [get(m,'yLabels'),get(m,'xLabels')];
le = legend(labels{:});
set(le,'interpreter','none');
grfun.movetosubplot(le,1,3,3);

%% Help on IRIS Functions Used in This File
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in an HTML browser window.
%
%    help model/subsasgn
%    help model/chksstate
%    help model/solve
%    help model/sstatedb
%    help model/simulate
%    help model/lhsmrhs
%    help plan
%    help plan/plan
%    help plan/nonlinearise
%    help dbase/dbextend
%    help dbase/dbplot
%    help grfun/ftitle
%    help grfun/bottomlegend
%    help grfun/movetosubplot
